
* dependencies:
*   estout
*   reghdfe and julia
*   blindschemes
*   grc1leg2
* churn


cap cd "/Users/Amit/OneDrive - Center for Disease Dynamics, Economics & Policy/RD response/dec dataverse_files/"

graph set window fontface "LM Roman 12"
set scheme plotplain // plottig

  use "Analysis data", clear

  recode General_Education* headedu* (1/7 = 1)  // http://microdata.gov.in/nada43/index.php/catalog/127/datafile/F4/V161
  recode Relation_to_Head* (5 = 3) (. 4 8 9 = 0)  // http://microdata.gov.in/nada43/index.php/catalog/127/datafile/F4/V157
  gen byte married = 2.Marital_Status

  bysort HHID (PID): gen headedu    = General_Education[1]
  by     HHID      : gen headage    = Age[1]
  by     HHID      : gen femalehead = !Male[1]
  recode headedu (. = 0)

  preserve
use "/Users/Amit/OneDrive - Center for Disease Dynamics, Economics & Policy/RD response/2024/Analysis data/nfhs4_migrant_clean.dta", clear
  probit migrant Age Agesq nevermarried child grandchild mpceq2 mpceq3 mpceq4 mpceq5 female head wifehead headage femalehead rural sc st obc muslim christian sikh HH_Size if work & birthyear>=1985 & birthyear<=1990
  restore
  gen double migrant85_90SNB = normprob(            _b[female]*femaleSNB +_b[Age]*AgeSNB +_b[Agesq]*AgeSNB^2	+_b[nevermarried]*nevermarriedSNB	+_b[head]*headSNB	+_b[wifehead]*wifeheadSNB	+_b[child]*childSNB	+_b[grandchild]*grandchildSNB	+_b[headage]*headageSNB	+_b[femalehead]*femaleheadSNB	+_b[rural]*ruralSNB	+_b[sc]*scSNB	+_b[st]*stSNB	+_b[obc]*obcSNB	+_b[muslim]*muslimSNB	+_b[christian]*christianSNB	+_b[sikh]*sikhSNB	+_b[HH_Size]*HH_SizeSNB	+_b[mpceq2]*mpceq2SNB	+_b[mpceq3]*mpceq3SNB	+_b[mpceq4]*mpceq4SNB	+_b[mpceq5]*mpceq5SNB)
  gen double migrant85_90    = normprob(_b[_cons] + _b[female]*femaleSNB +_b[Age]*AgeSNB +_b[Agesq]*AgeSNB^2	+_b[nevermarried]*nevermarriedSNB	+_b[head]*headSNB	+_b[wifehead]*wifeheadSNB	+_b[child]*childSNB	+_b[grandchild]*grandchildSNB	+_b[headage]*headageSNB	+_b[femalehead]*femaleheadSNB	+_b[rural]*ruralSNB	+_b[sc]*scSNB	+_b[st]*stSNB	+_b[obc]*obcSNB	+_b[muslim]*muslimSNB	+_b[christian]*christianSNB	+_b[sikh]*sikhSNB	+_b[HH_Size]*HH_SizeSNB	+_b[mpceq2]*mpceq2SNB	+_b[mpceq3]*mpceq3SNB	+_b[mpceq4]*mpceq4SNB	+_b[mpceq5]*mpceq5SNB)

  replace Date_of_Survey = date(string(Date_of_Survey, "%06.0f"), "DM20Y")
  replace Date_of_Survey = Date_of_Survey - 366 if Date_of_Survey > td(30jun2012)  // Amit Summan suggestion: pull a few post-survey observations back a year
  format %td Date_of_Survey
  gen int Month_of_Survey = mofd(Date_of_Survey)

  gen svyyear = yofd(Date_of_Survey)
  gen int birthyear    = svyyear    - Age
  gen int birthyearSNB = svyyearSNB - AgeSNB
  gen byte TSNB = birthyearSNB >= uipSNB
  replace uip = round(uip, 1)
  recast int uip
  replace uip = . if uip_SD > .1
  gen byte T = birthyear >= uip

  gen double lwage = ln(Wages) if !inlist(Usual_Principal_Activity_Status, 11,12,21)  // don't count self-employment wages
  gen double ldwage = lwage - ln(Days_worked)  // convert to wages per time worked

 /* compute CPIs as function of dates, interpolating from monthly data, 6/2011-12/2012, from https://data.imf.org/?sk=fffaba71-e02b-4e2c-87bc-45cf050eb5bc&hide_uv=1
AS: additional measure of inflation added here, suffixed with WB, HCPI from
https://www.worldbank.org/en/research/brief/inflation-database
  */
  * Step 1: Generate CPI variables using IMF, World Bank (WB)
gen CPI = .                    // IMF CPI
gen CPI_WB = .                 // World Bank CPI
gen byte hasdate = svyyear < .

mata {
    // IMF CPI Data
    CPI = 107.44 \ 109.71 \ 110.28 \ 111.98 \ 112.55 \ 113.12 \ 112.63 \ 113.02 \ 113.62 \ 114.51 \ 
          116.10 \ 117.19 \ 118.57 \ 120.36 \ 121.85 \ 122.94 \ 123.83 \ 124.32 \ 124.52

    // World Bank CPI Data
    CPI_WB = 72.30 \ 73.83 \ 74.21 \ 75.36 \ 75.74 \ 76.12 \ 75.36 \ 75.74 \ 76.12 \ 76.89 \ 
             78.42 \ 78.80 \ 79.57 \ 81.10 \ 81.86 \ 82.24 \ 83.01 \ 83.39 \ 83.77



    // Compute differences for interpolation
    dCPI = CPI[|2\.|] - CPI[|.\length(CPI)-1|]
    dCPI_WB = CPI_WB[|2\.|] - CPI_WB[|.\length(CPI_WB)-1|]

    // Link to survey data
    st_view(surveyCPI=., ., "CPI", "hasdate")
    st_view(surveyCPI_WB=., ., "CPI_WB", "hasdate")
    date = st_data(., "Date_of_Survey", "hasdate")
    
    y = year(date); m = month(date); d = day(date); my = mofd(date)
    i = y * 12 + m :- (2011*12+6)  // indexes in CPI tables for last value before date

    // IMF CPI interpolation
    surveyCPI[,] = CPI[i] + d:/(dofm(my:+1) - dofm(my)) :* dCPI[i]

    // World Bank CPI interpolation
    surveyCPI_WB[,] = CPI_WB[i] + d:/(dofm(my:+1) - dofm(my)) :* dCPI_WB[i]


}

* Step 2: Adjust variables using IMF CPI
gen double lwagereal = lwage - ln(CPI/100)


* Step 3: Adjust variables using World Bank CPI
gen double lwageWBreal = lwage - ln(CPI_WB/100)


  sum WEIGHT, detail
  gen WTtr = min(WEIGHT, r(p50)+5*(r(p75)-r(p50)))  // clip extreme weights to median + 5 * IQR (Potter and Zheng 2015)

save analysis_datav2,replace

* best replication; then monthly survey effects instead; then inflation adjusted
cap erase esttab.rtf
 egen fe_id=group (Age districtcode_nsso)

global controlsSNB urban Male marriedSNB HH_Size headage i.Social_Group i(2/4)bn.Religion i(1 2 3 6)bn.Relation_to_Head femalehead i(8 10 12 13)bn.(General_Education headedu) migrant85_90SNB
global controls    urban Male married    HH_Size headage i.Social_Group i(2/4)bn.Religion i(1 2 3 6)bn.Relation_to_Head femalehead i(8 10 12 13)bn.(General_Education headedu) migrant85_90
global controlsSNB_noed urban Male marriedSNB HH_Size headage i.Social_Group i(2/4)bn.Religion i(1 2 3 6)bn.Relation_to_Head femalehead i(8 10 12 13)bn.( headedu) migrant85_90SNB

***Table 1 results replicate DR main results and then focus on related specifications 
{
eststo clear
    eststo DR_SNB: reghdfe lwagereal TSNB $controlsSNB if inrange(birthyearSNB,1985,1990) [aw=WTtr], cluster(District_Code) a(Age#District_Code) keepsing
	eststo DR_SNB_WB: reghdfe lwageWBreal TSNB $controlsSNB if inrange(birthyearSNB,1985,1990) [aw=WTtr], cluster(District_Code) a(Age#District_Code) keepsing
	  eststo SNB_m: xtreg lwage TSNB $controlsSNB  if inrange(birthyearSNB,1985,1990) & Male==1, fe i(fe_id) robust cluster(districtcode_nsso) 
	   eststo DR_SNB_m: reghdfe lwagereal TSNB $controlsSNB if inrange(birthyearSNB,1985,1990) & Male==1 [aw=WTtr], cluster(District_Code) a(Age#District_Code) 
	    eststo DR_SNB_m_WB: reghdfe lwageWBreal TSNB $controlsSNB if inrange(birthyearSNB,1985,1990) & Male==1 [aw=WTtr], cluster(District_Code) a(Age#District_Code) keepsing
		eststo DR_SNB_m_WB_old: reghdfe lwageWBreal TSNB $controlsSNB if inrange(birthyearSNB,1980,1990) & Male==1 [aw=WTtr], cluster(District_Code) a(Age#District_Code) keepsing

	 
**AS: formatting of table, adding significance indicators, coverting to percentage change
// List of all models
local models DR_SNB DR_SNB_WB SNB_m DR_SNB_m DR_SNB_m_WB DR_SNB_m_WB_old 
*DR DR_WB SNB_m2 DR_m DR_m_WB DR_m_WB_old

// Create a numeric matrix to store results: Coefficient, P-value, and Observations
matrix results = J(3, `:word count `models'', .)
local rownames "Coefficient P_Value Observations"

// Loop through models to extract and transform results
local col = 1
foreach model in `models' {
    est restore `model'

    // Determine the variable of interest (TSNB or T)
    local var_interest = cond(strpos("`:colnames e(b)'", "TSNB"), "TSNB", "T")
    
    // Extract the coefficient and apply the transformation
    local coef = exp(_b[`var_interest']) - 1

    // Round the coefficient to two decimal places
    local rounded_coef = round(`coef', 0.001)

    // Calculate the p-value using the t-distribution
    local se = sqrt(e(V)[colnumb(e(b), "`var_interest'"), colnumb(e(b), "`var_interest'")])
    local tstat = _b[`var_interest'] / `se'
    local df = e(df_r)  // Residual degrees of freedom
    local pvalue = 2 * ttail(`df', abs(`tstat'))

    // Extract number of observations
    local obs = e(N)

    // Store numeric values into the matrix
    matrix results[1, `col'] = `rounded_coef' // Transformed Coefficient
    matrix results[2, `col'] = `pvalue'       // P-Value
    matrix results[3, `col'] = `obs'          // Observations

    local col = `col' + 1
}

// Add row and column names to the matrix
matrix rownames results = Coefficient P_Value Observations
matrix colnames results = `models'

// Safely close the file handle if it already exists
capture file close output

// Open the file for manual export
file open output using "results.csv", write replace

// Write the header row
file write output "Statistic"
foreach model in `models' {
    file write output ",`model'"
}
file write output _n

// Write the results row by row
forval i = 1/3 {
    local rowname: word `i' of `rownames'
    file write output "`rowname'"
    forval j = 1/`=colsof(results)' {

        if `i' == 1 {
            // coefficient: 3 decimals + stars
            local val  : display %9.3f results[`i',`j']
            local pval = results[2,`j']
            local stars ""
            if      `pval' < 0.01 local stars = "**"
            else if `pval' < 0.05 local stars = "*"
            else if `pval' < 0.10 local stars = "+"

            file write output ",`val'`stars'"
        }
        else if `i' == 2 {
            // p-value: 3 decimals
            local val : display %9.3f results[`i',`j']
            file write output ",`val'"
        }
        else {
            // observations: integer (with commas); use %9.0f if you don't want commas
local val : display %9.0f results[`i',`j']
            file write output ",`val'"
        }
    }
    file write output _n
}


// Close the file
file close output
}

***Table A2 results replicate DR appendix results and then focus on related specifications 
{
 eststo DR: reghdfe lwagereal T $controls if inrange(birthyear,1985,1990) [aw=WTtr], cluster(District_Code) a(Age#District_Code) 
eststo DR_WB: reghdfe lwageWBreal T $controls if inrange(birthyear,1985,1990) [aw=WTtr], cluster(District_Code) a(Age#District_Code) 
eststo SNB_m: xtreg lwage TSNB $controlsSNB  if inrange(birthyearSNB,1985,1990) & Male==1, fe i(fe_id) robust cluster(districtcode_nsso) 
eststo DR_m: reghdfe lwagereal T $controls if inrange(birthyear,1985,1990) & Male==1 [aw=WTtr], cluster(District_Code) a(Age#District_Code) 
eststo DR_m_WB: reghdfe lwageWBreal T $controls if inrange(birthyear,1985,1990) & Male==1 [aw=WTtr], cluster(District_Code) a(Age#District_Code) 
eststo DR_m_WB_old: reghdfe lwageWBreal T $controls if inrange(birthyear,1980,1990) & Male==1 [aw=WTtr], cluster(District_Code) a(Age#District_Code) 


**AS: formatting of table, adding significance indicators, coverting to percentage change
// List of all models
local models DR DR_WB SNB_m DR_m DR_m_WB DR_m_WB_old

// Create a numeric matrix to store results: Coefficient, P-value, and Observations
matrix results = J(3, `:word count `models'', .)
local rownames "Coefficient P_Value Observations"

// Loop through models to extract and transform results
local col = 1
foreach model in `models' {
    est restore `model'

    // Determine the variable of interest (TSNB or T)
    local var_interest = cond(strpos("`:colnames e(b)'", "TSNB"), "TSNB", "T")
    
    // Extract the coefficient and apply the transformation
    local coef = exp(_b[`var_interest']) - 1

    // Round the coefficient to two decimal places
    local rounded_coef = round(`coef', 0.001)

    // Calculate the p-value using the t-distribution
    local se = sqrt(e(V)[colnumb(e(b), "`var_interest'"), colnumb(e(b), "`var_interest'")])
    local tstat = _b[`var_interest'] / `se'
    local df = e(df_r)  // Residual degrees of freedom
    local pvalue = 2 * ttail(`df', abs(`tstat'))

    // Extract number of observations
    local obs = e(N)

    // Store numeric values into the matrix
    matrix results[1, `col'] = `rounded_coef' // Transformed Coefficient
    matrix results[2, `col'] = `pvalue'       // P-Value
    matrix results[3, `col'] = `obs'          // Observations

    local col = `col' + 1
}

// Add row and column names to the matrix
matrix rownames results = Coefficient P_Value Observations
matrix colnames results = `models'

// Safely close the file handle if it already exists
capture file close output

// Open the file for manual export
file open output using "results_appendix.csv", write replace

// Write the header row
file write output "Statistic"
foreach model in `models' {
    file write output ",`model'"
}
file write output _n

// Write the results row by row
local rownames "Coefficient P_Value Observations"
forval i = 1/3 {
    local rowname: word `i' of `rownames'
    file write output "`rowname'"
    forval j = 1/`=colsof(results)' {
        if `i' == 1 {
            // coefficient: 3 decimals + stars
            local val  : display %9.3f results[`i',`j']
            local pval = results[2,`j']
            local stars ""
            if      `pval' < 0.01 local stars = "**"
            else if `pval' < 0.05 local stars = "*"
            else if `pval' < 0.10 local stars = "+"
            file write output ",`val'`stars'"
        }
        else if `i' == 2 {
            // p-value: 3 decimals
            local val : display %9.3f results[`i',`j']
            file write output ",`val'"
        }
        else {
            // observations: integer (no thousands separator to keep CSV clean)
            local val : display %9.0f results[`i',`j']
            file write output ",`val'"
        }
    }
    file write output _n
}


// Close the file
file close output
}



***
*** figure 1: change in treatment obs ove survey period
***

use analysis_datav2, clear

// Step 1: Adjust survey months to start at July 2011
gen Survey_Calendar_Month = Month_of_Survey - 617  // 618 becomes 1 (July 2011)

// Step 2: Collapse data with weights and filter by birth years
collapse (mean) TSNB [aw=WTtr] if birthyearSNB >= 1985 & birthyear <= 1990, by(Survey_Calendar_Month) 


// Generate the graph with the correlation coefficient as a note
twoway (scatter T Survey_Calendar_Month, mcolor(blue)) /// Scatter plot
       (lfit T Survey_Calendar_Month, lcolor(red)), /// Fitted line
       ytitle("Proportion Treated", margin(medium)) /// Increase gap between y-axis and label
       xtitle("") /// Remove x-axis title
       ylab(0(0.1).8, angle(0)) /// Adjust y-axis labels
       xlabel(1 "07-11" 2 "08-11" 3 "09-11" ///
              4 "10-11" 5 "11-11" 6 "12-11" ///
              7 "01-12" 8 "02-12" 9 "03-12" ///
              10 "04-12" 11 "05-12" 12 "06-12", ///
              angle(0) labsize(small)) /// Rotate x-axis labels and remove ticks
       legend(off) ///
       graphregion(color(white))

// Export the graph as a PNG file
graph export "proportion_treated.png", as(png) replace width(2000) height(1500)






/**************************************************************************
 RI do-file (DR model 3 used in randomization exercise)
 Produces: ri.xlsx with sheets main_sim1..4
   Each sheet has 3 columns (b_sim, t, ri) for `d' draws.
**************************************************************************/


* ----- settings -----
local outxls "/Users/Amit/OneDrive - Center for Disease Dynamics, Economics & Policy/RD response/dec dataverse_files/ri.xlsx"
local d_sim = 1000   // churn length to build sim columns
local d     = 1000    // draws used in regressions (rows r1..r`d')

* Controls (define once)
global controlsSNB ///
    urban Male marriedSNB HH_Size headage ///
    i.Social_Group i(2/4)bn.Religion i(1 2 3 6)bn.Relation_to_Head ///
    femalehead i(8 10 12 13)bn.(General_Education headedu) ///
    migrant85_90SNB

/**************************************************************************
 Build simulation files rT1, rT2, rT3, rT4
 - rT1: churn uipSNBsim across districts (later convert to TSNB via threshold) keep distibution across years same
 - rT2: churn TSNB within each birthyear, then append across 1984..1991
 - rT3: churn TSNB for birthyears 1986..1989
 - rT4: churn TSNB for birthyears 1985..1990
**************************************************************************/
set seed 123456789

* rT1
use analaysis_data2, clear
gen uipSNBsim = uipSNB
duplicates drop districtcode_nssoSNB, force
keep *district* uip*
drop if missing(uipSNBsim)
forvalues x = 1/`d_sim' {
    churn uipSNBsim, clear
    gen uipSNBsim`x' = uipSNBsim
}
save rT1.dta, replace

* rT3
use analaysis_data2, clear
gen byte TSNBsim = TSNB
duplicates drop districtcode_nssoSNB birthyearSNB, force
keep *district* TSNB uip* birthyearSNB
drop if missing(TSNB)
forvalues x = 1/`d_sim' {
    churn TSNB if inrange(birthyearSNB,1986,1989), clear
    gen TSNBsim`x' = TSNB
}
save rT3.dta, replace

* rT4
use analaysis_data2, clear
gen byte TSNBsim = TSNB
duplicates drop districtcode_nssoSNB birthyearSNB, force
keep *district* TSNB uip* birthyearSNB
drop if missing(TSNB)
forvalues x = 1/`d_sim' {
    churn TSNB if inrange(birthyearSNB,1985,1990), clear
    gen TSNBsim`x' = TSNB
}
save rT4.dta, replace

* rT2 (equal counts within each birthyear; churn within birthyear, then append)
forvalues y = 1984/1991 {
    preserve
    use analaysis_data2, clear
    duplicates drop districtcode_nssoSNB birthyearSNB, force
    keep if birthyearSNB==`y'
    keep *district* TSNB uip* birthyearSNB
    drop if missing(TSNB)
    forvalues x = 1/`d_sim' {
        churn TSNB, clear
        gen TSNBsim`x' = TSNB
    }
    save rTv2`y'.dta, replace
    restore
}
use rTv21984.dta, clear
forvalues y=1985/1991 {
    append using rTv2`y'.dta
}
save rT2.dta, replace

/**************************************************************************
 Run regressions for z=2,3,4 (merge TSNBsim columns from rTz),
 and z=1 (build TSNBsim from uipSNBsim threshold).
 Write ONLY the 3 needed columns (coef, p, ri) to ri.xlsx:
   main_sim`z'
**************************************************************************/

* ---------- z = 2/4 ----------
forvalues z = 2/4 {
    use analaysis_data2, clear
    gen byte TSNBsim = TSNB
    merge m:1 districtcode_nssoSNB birthyearSNB using rT`z'.dta
    keep if _merge==3
    drop _merge

    * Preallocate 3-column matrix (b_sim, p, ri) with fixed names
    mat main_sim`z' = J(`d', 3, .)
    local rn
    forvalues i=1/`d' {
        local rn `rn' r`i'
    }
    matrix colnames main_sim`z' = c1 c2 c3
    matrix rownames main_sim`z' = `rn'

    forvalues x = 1/`d' {
        * Share of unchanged treatment vs sim
        cap drop T_diff
        gen T_diff = TSNB == TSNBsim`x'
        replace T_diff = 0 if missing(T_diff)
		 quietly  sum T_diff [aw=WTtr] if inrange(birthyearSNB,1985,1990)
        local share = r(mean)

        * RD model 3 
        quietly reghdfe lwage TSNBsim`x' $controlsSNB ///
            if inrange(birthyearSNB,1985,1990) [aw=WTtr], ///
            cluster(District_Code) a(Age#District_Code) keepsing
        local b = _b[TSNBsim`x']
        local p = 2*ttail(e(df_r), abs(_b[TSNBsim`x']/_se[TSNBsim`x']))

        mat main_sim`z'[`x',1] = `b'
        mat main_sim`z'[`x',2] = `p'
        mat main_sim`z'[`x',3] = `share'
    }
	
	

    * Write exactly three columns; keep row/col names for layout
    putexcel set "`outxls'", sheet("main_sim`z'", replace) modify
	*nformat(number_d2)
    putexcel A1 = matrix(main_sim`z'), names 
}

* ---------- z = 1 (TSNBsim from uipSNBsim threshold) ----------
use analaysis_data2, clear
gen uipSNBsim = uipSNB
merge m:1 districtcode_nssoSNB using rT1.dta
keep if _merge==3
drop _merge

mat main_sim1 = J(`d', 3, .)
local rn
forvalues i=1/`d' {
    local rn `rn' r`i'
}
matrix colnames main_sim1 = c1 c2 c3
matrix rownames main_sim1 = `rn'

forvalues x = 1/`d' {
    cap drop TSNBsim`x'
    gen TSNBsim`x' = birthyearSNB >= uipSNBsim`x'

    cap drop T_diff
    gen T_diff = TSNB == TSNBsim`x'
    replace T_diff = 0 if missing(T_diff)
	quietly sum T_diff [aw=WTtr] if inrange(birthyearSNB,1985,1990)
    local share = r(mean)
	
quietly reghdfe lwage TSNBsim`x' $controlsSNB ///
        if inrange(birthyearSNB,1985,1990) [aw=WTtr], ///
        cluster(District_Code) a(Age#District_Code) keepsing
    local b = _b[TSNBsim`x']
    local p = 2*ttail(e(df_r), abs(_b[TSNBsim`x']/_se[TSNBsim`x']))
	


    mat main_sim1[`x',1] = `b'
    mat main_sim1[`x',2] = `p'
    mat main_sim1[`x',3] = `share'
}

putexcel set "`outxls'", sheet("main_sim1", replace) modify
putexcel A1 = matrix(main_sim1), names nformat(number_d2)


/**************************************************************************
 test and graph 
**************************************************************************/

clear
set more off



local d = 1000
local end = `= `d' + 1'   // +1 for the header row

forvalues z = 1/4 {
    // 1) Import and prep
    import excel ///
        "/Users/Amit/OneDrive - Center for Disease Dynamics, Economics & Policy/RD response/dec dataverse_files/ri.xlsx", ///
        sheet("main_sim`z'") ///
        firstrow clear cellrange(A1:D`end')

    rename (c1 c2 c3) (b_sim t ri)
    drop if missing(b_sim)

    // 2) Observed coefficient
    scalar b_obs = 0.119
    local point_estimate = b_obs

    // 3) Mean of simulated coefficients
    egen mean_sim = mean(b_sim)
    local randomization_mean = mean_sim

    // 4) Two-tailed p-value
    gen farther = abs(b_sim - mean_sim) >= abs(b_obs - mean_sim)
    summarize farther
    scalar p_twotail = r(sum) / _N

    // 5) % of unchanged treatment obs (rounded display)
    summarize ri
local change_display = string(100 * r(mean), "%4.2f")

    // 6) Plot
    twoway (kdensity b_sim), ///
        xline(`point_estimate',     lcolor(red)   lpattern(solid)) ///
        xline(`randomization_mean', lcolor(black) lpattern(dash))  ///
     xlabel( -0.20 "-0.20"  -0.10 "-0.10" 0 "0" 0.10 "0.10"  0.20 "0.20"  ///
               `point_estimate'     `" " ""PE" "' ///
               `randomization_mean' `" " ""RM" "' ) ///
        xscale(r(-0.25 0.25)) ///
		graphregion(color(white)) ///
        yscale(off) ///
        xtitle("") ///
        title("       Model `z'", size(medium) justification(left)) ///
        note("`change_display'% of treatment obs unchanged" "p = `:display %4.3f p_twotail'", size(small) position(13)) ///
        name(graph`z', replace)

    graph save "graph`z'.gph", replace
    graph export "graph`z'.png", replace
}

* Combine the four  graphs
graph combine "graph1.gph" "graph2.gph" "graph3.gph" "graph4.gph", ///
    graphregion(color(white)) plotregion(color(white))

graph save "combined_graph_with_legend.gph", replace
graph export "combined_graph_with_legend.png" , replace width(3000)
